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THEORY  AND  COMPUTATION  OF  MULTILAYER 
COMPOSITES 


E.  STEIN  AND  J.  TESSMER 
Institute  for  Structural  Mechanics  and 
Computational  Mechanics,  University  of  Hannover, 
Appelstrafie  9A,  D-S0167  Hannover,  Germany 


Abstract.  An  hierarchical  concept  for  the  analysis  of  thin-walled  com- 
posite shells  is  presented  within  the  finite  element  method  for  nonlinear 
deformation  of  thin-walled  composite  structures.  For  non-disturbed  subdo- 
mains of  the  structure  an  effective  4-node  shell  element  with  5 or  6 d.o.f. 
per  node  is  used  for  the  whole  laminate.  For  the  analysis  of  3D-stress  states 
a multidirector  shell  element  along  the  thickness  with  piecewise  polynomi- 
als is  presented,  n physical  layers  are  approximated  by  N numerical  layers, 
discretized  with  hierarchical  trial-  and  test  functions.  The  coupling  of  both 
elements  is  performed  by  special  transition  elements.  An  example  shows 
the  technique  and  efficiency  of  coupling  both  types  of  elements. 


1.  Introduction 

Composites  are  typically  used  for  light-weight  structures,  and  in  many 
engineering  fields  there  are  attempts  to  replace  components  with  clas- 
sical materials  (steel,  concrete)  by  fiber  reinforced  materials.  Therefore 
we  focus  on  the  computation  of  thin-walled  laminated  composites.  Due 
to  anisotropy  and  inhomogenity  these  structures  show  rather  complicated 
states  of  stresses  and  strains  such  that  several  adequate  mechanical  models 
for  thin  composite  structures  have  been  developed  in  recent  years,  see  [5,  9]. 

Different  F E-methods  were  applied  with  respect  to  special  purposes 
(global  deformation  analysis,  stress  analysis,  first  ply  failure  criteria,  local 
stress  singularities,  crack  analysis),  see  e.g.  [2,  8].  By  admitting  a priori 
delaminations  between  plies  and  regarding  their  growth,  load  depending 
disturbances  of  the  perfect  geometry  of  cross-sections  are  growing  and  di- 
minish the  critical  load  within  stability  analysis.  We  present  an  hierarchical 
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multi-layer  shell  model  as  a framework  for  the  efficient  nonlinear  analysis 
of  delamination  and  failure  process,  focusing  on  basic  topics  of  composite 
theory  and  computation.  In  [10]  the  failure  analysis  for  this  model  is  shown. 
For  calculating  a complete  3D  stress  state  in  disturbed  subdomains  a mul- 
tidirector shell  kinematic  with  deformation  modes  in  thickness  direction  is 
used  and  implemented  in  FEM.  Within  this  kinematics,  the  interpolation 
in  thickness  direction  with  hierarchical  polynomials  per  layer  is  indepen- 
dent from  in  plane  interpolation.  Since  dominant  parts  of  composite  shells 
usually  reside  in  2D  stress  states  with  little  transvere  shear  deformation,  a 
conventional  Reissner-Mindlin  kinematic  is  used  for  regular  parts.  In  these 
areas  normal  stresses  S33  are  much  smaller  than  normal  in-plane  stresses. 
Therefore,  S33  is  neglected  in  the  classical  laminate  theory.  For  coupling 
both  kinematic  types  within  a finite  element  method  a transition  element 
is  applied. 

2.  Kinematics 

The  considered  thin-walled  composite  structures  consist  of  a layerwise  build 
up  with  n physical  layers  j of  thickness  h3  and  N numerical  layers  i which 
can  collect  or  subdivide  the  physical  layers  for  numerical  calculations,  fig. 
1 and  2 

The  position  vector  Xo  of  the  reference  surface  So  is  parametrizised  by 
convected  coordinates  0a.  An  orthonormal  basis  system  tfc(0Q)  is  attached 
to  this  surface  where  t3  is  a normal  vector  and  03  the  coordinate  in  thick- 
ness direction.  The  transformation  between  the  different  base  systems  is 
given  by  tfc(0")  = Ro(0a)  e*,  where  Ro  is  a proper  orthogonal  tensor. 

Depending  on  the  desired  accuracy  and  numerical  effort  two  different 
types  of  kinematics  are  applied.  One  is  the  standard  shell  kinematic  with 
one  director,  yielding  a 2D-FE-formulation  with  5 or  6 d.o.f.  per  node, 
and  the  other  is  a 3D-multi-director  kinematic  which  still  yields  a 2D-like 
finite  element  data  structure  by  a priori  C°-continuous  hierarchical  shape 
functions  in  thickness  direction,  such  we  get  the  numerical  complexity  of  a 
3D-brick  formulation. 

2.1.  INEXTENSIBLE  ONE-DIRECTOR  KINEMATICS 

The  position  vectors  of  the  reference  and  the  current  configuration  of  the 
shell  body  are  given  by 

X(0“,  ©3)  = Xo(©a)  + ©3t3(0Q),  -hu<e3<h0,  m 

x(©“,©3)  = xo(©a)  + 03d(0“),  W 
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with  xo(0“)  = Xo(0“)  + uo(0Q) , 

d = a3  = Rt3,  ||dji  = ||t3||  = 1,  (9) 

R = /(A/3)  £ SO 3 , A/3  = rotational  vector , ' ' 

v = [u0;A (3\T . 

uo  is  the  displacement  vector  of  the  reference  surface,  d the  director  vec- 
tor, and  R is  an  orthogonal  tensor  generated  by  application  of  the  Euler- 
Rodriguez-formula  on  the  rotational  vector  A/3,  see  e.g.  [1]. 

The  Green-Lagrangian  strain  tensor  E = EijG1  <g>  G-7  is  formulated 
with  respect  to  the  tangential  base  vectors  Gj,  = X,*,  = dX/dOk  for  the 
inextensible  thinwalled  structure. 

Eij  = |(x,j-x,j  -X,jX,j ) = £ij  + k,j©3  + | ...(03)2 , (3) 

with  the  approximations 

£aP  — 0.5(x,a  X,^3  X,a  X,f, ) , 

7a  = 2sa3  = (xo,a  ' d Xo ,a  ' ^3)  ? (4) 

— 0-b[(xo,a  ’ 'Xo,^)  (Xo^q  ' ^3)/?  ~ 1“ ^3 ,q.  -Xo,/?)]  • 

2.2.  EXTENSIBLE  MULTIDIRECTOR  KINEMATICS 

For  the  description  of  3D-stress-states  an  extended  multidirector  kinematic 
is  used.  The  position  vectors  of  the  reference  and  the  current  configuration 
of  the  shell  body  are  given  by 

X(0“,  ©3)  = Xo(©“)  + 03t3(0a)  -hu<03<h0, 

x(©“,03)  = X(©“,©3)  + u(0“,  03) . U 
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A multiplicative  decomposition  of  u is  introduced  as 

u(0a,03)  = £ u;(0a)  wi(Q3)  = #(©3)  u(©“) , 

1=1 

u(0“)  = [uX,  U2,  ...  , Um,  ...  , Us+1,  Us+2,  ...  , U*+Tn,  ...  , UN(m_1)+1]T  (6) 
a1  (layer  X)  n>  (layer  »') 

where  k is  the  number  of  interpolation  points  over  the  thickness,  wi  are 
shape  function  (piecewise  linear  for  k = N -f  1),  s — (i  — 1)  (m  — 1)  is  a 
layer  index,  and  m is  the  polynomial  degree,  (2  < m < 4).  In  thickness 
direction  hierarchical  shape  functions  are  defined  for  each  numerical  layer 
i as  #'(03)  . 

ul(0a,  03)  = 4>‘(©3)  ul(0“) , u*(©Q)  = [uj,  u2>  ...  , u^]T , 

«i  = [»i„  «i„,  u[z]T , = [4>\i,  41,  ...  , 4L i]  , 

4(C)  = 0.5(1  -C)>  4(C)  = 1-C2,  C = 203At- 

4(C)  = 0.5(1  + C) , 4(C)  = (1  - C2)C  ■ 

This  leads  to  an  a priori  discretization  in  thickness  direction  within  the 
analytical  shell  theory,  analogous  to  the  analytical  thickness  integration  of 
stress  resultants  in  classical  shell  models.  These  shape  functions  are  inde- 
pendent from  the  in-plane  trial-functions,  see  sec.  5.2.  For  linear  plates  see 
[8,  9].  A layerwise  linear  interpolation  is  realized  by  using  cf>\  and  <j>\  only. 

In  contrast  to  eq.(4)  the  Green-Lagrangian  strain  tensor  is  computed 
and  applied  in  its  complete  form  , eq.(3). 
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3.  Material  Equations  on  Macro  Level 

3D-material  equations  are  formulated  for  each  physical  layer  of  the  com- 
posite laminate  . These  have  to  be  transformed  from  the  fiberoriented  basis 
to  the  local  basis  of  the  shell,  see  [5].  In  case  of  standard  shell  kinematics 
they  are  integrated  a priori  over  the  thickness  of  the  overall  thickness. 


3.1.  UNI-DIRECTIONAL  TRANSVERSAL  ISOTROPIC  LAYERS 


The  constitutive  equations  for  unidirectional  (UD)  layers  of  laminated  struc- 
tures are  derived  under  the  assumption  of  small  strains  and  the  validity  of 
St.  Venant  Kirchhoff  material  equations  with  linearized  strains. 

Hence,  the  free  energy  TV(E)  = |E  • CE  is  formulated  as  a quadratic 
function  of  E.  The  material  tensor  C is  constant  with  respect  to  the  in- 
variants trE,  tr(E2),  ao  • Eao,  ao  • E2ao.  Using  the  elastic  parameters 
A,  fiTi  //£,,  <*,/?,  the  free  energy  W,  the  stresses  S and  the  constant  material 
tensor  C follow  as 


W(E) 


with 


|A(trE)2  + //Ttr(E2)  + a(ao  ■ Eao)trE 
+2 (fiL  ~ /-tr)(ao  • E2ao)  + |/3(a0  • Ea0)2, 


S = *fP  = C.E, 


(8) 


For  UD-layers  with  the  fiber  direction  ao  = ei  = (1,0,0)  and  by  compu- 
tation of  CC_1  = 1 the  components  of  the  compliance  matrix  €_1  result 
with  the  common  elasticity  constants  in 


cr1 


1 /E\  ~~vi2/E\  ~vi2/Ei  0 0 0 

—V12/E1  1/  E2  —V23/E2  0 0 0 

~vi  2/ Ei  —V23/E2  I/E2  0 0 0 

0 0 0 1/G12  0 0 

0 0 0 0 I/G12  0 

0 0 0 0 0 I/G23 


3.2.  HYPERELASTIC  ISOTROPIC  INTERMEDIATE  LAYERS 

For  composite  structures  with  macroscopicly  seperated  parts  of  anisotropic 
and  isotropic  material,  e.g.  tires,  also  a hyperelastic  isotropic  material  equa- 
tion is  formulated  with  respect  to  the  invariants  of  the  Cauchy-Green  tensor 
C = FtF  = 2E  + 1.  The  stored  energy  by  ‘Blatz  & Ko’,  [7],  yields 

W = — 3)  + (1  — /)(J2  - 3)+ 

^ 1 _ 9,,  (10) 

t1— ty/MT**  - 1)  + (i  - /)«,■-’•’  - 1)». 

Vq 
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(11) 


With  j\  = h , j2  = h/h,  J3  = Vh  = J\ 


h = trC  , h - 0.5[(trC)2  - trC2]  , J3  = detC  . 
The  2.  Piola-stress  tensor  S and  the  material  tensor  C follow  as 
dW  d 

S = ~8E  =^QcW(h(C),I2(C),I3(C)), 

C = 4^W(Il(C),I2(C),I3(C)). 

=>  s = /xotz+^a-zm-^ia-^c+^c-1, 

C = 2^(1-/)(1®1-I)  + 2Mo/?i^ 

+ [Mo  ^ ~ 2u0^2  ^°^"(1  ~ /)]^ 

+ 2j~(l  - /)(C“l  ® C + C ® C"1) 

- 27^(1  - /)(C-X  ®1  + 1®  C"1)  , 

Zt'n 

with  /3j  = + (1  - /)J  l-2-o  - J2(l  - /)  , 


02 

dC-1 

dC 


— 2kq  Zt'o 

/Jl-2-0  +(l-/)Jl-2^0  , 

-C^ICT1. 


(12) 


3.3.  RESULTANT  FORMULATION  FOR  ONE-DIRECTOR  SHELL 


Since  the  material  equations  for  one-director  shells  are  formulated  in  resul- 
tants of  stresses  and  strains,  an  adequate  transformation  and  integration 
of  the  material  law,  eq.(9),  is  conducted.  With  the  assumption  of  S33  = 0 
the  matrix  06x6  can  be  statically  kondensed  to  Cfx5  . Such,  we  get  for  the 
fiberorientated  basis 


Sj?  = Cp  Ep  with  c£  = 

cf  0 

0 Cf 

Cln  = Ej/(  1 ~ ^fe)  > 

Cml3  = 0 1 

Cfu  = On 

Cln  = E^/(  1-^f), 

C'm23  ~ 0 > 

c£2  = o, 

c£22  = E2/(  1 - ^f) . 

CZs3  = G 12, 

Cs22~  G23 

(13) 


SR=  [SU;S22;S12;S13;S23]  , 


ER  = [En;  E22‘i  2E12;  2Ei3;  2E23]  . 


For  the  integration  over  the  thickness  all  matrices  for  each  layer  j must  be 
transformed  to  the  local  basis  of  the  reference  surface,  see  [4]. 


Cf  = TjTCpTj 


c L o 
o Cj 


(14) 
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N Dm  D mb  0 £ 

=»  S = €-E=  M = Dlb  D6  0 • K , (15) 

. Q J o o dsJ  [7. 

with  S = [IV11;  IV22;  N12;  M11;  M22;  M12;  Q1;  £j?2]T  , 

E = [eu;  £22;  2ei2;  «n;  K22;  2k12;  71;  72]r 

Dm  = Dmi)  = 

D»  = ZjCitih  + gy,  d,  = EyCjW 


4.  Variational  Equations 

In  this  section  the  weak  form  of  equilibrium  (principle  of  virtual  work)  is 
given.  Linearization  yields  the  consistent  tangential  operator. 


4.1.  ONE-DIRECTOR  FORMULATION 


The  principle  of  virtual  work  is  given  in  the  material  description  for  the 
inextensible  shell  in  thickness  direction  and  external  loads  t at  the  reference 
surface,  using  (15) 

G(u,ri)  = JS-SEdV-  f t • 77 dlV  = 0 
(V)  (r<,) 

= f (Na/3S£Q0  + Mal3SKQp  + Qa6'ya)dQ  - Gext(i,r))  /16\ 

(ft)  1 ; 

= / S-SEdQ  - Gext(t,v)- 

(ft) 

For  constant  material  matrix  C and  conservative  loads  the  linearization 
(Gateaux-derivative)  of  (16)  yields 

DG(u,rj)  • Au  = / 6E-C-AEdQ  + f SASEdQ.  (17) 

(ft)  (ft)  ^ ’ 


4.2.  MULTIDIRECTOR  FORMULATION 

44 

The  laminate  structure  is  loaded  at  the  bottom-  and  top-surface  with 
p = pk  ek  . Hence,  the  virtual  work  reads  in  index  notation 

G{ u,  rj ) = f[  f Skl  SEki  JdO 3]  dQ2dQ1  - / p*  rjk  dtia  = 0 , n 
(ft)  (©3)  (ft.)  1 ; 

with  J(0‘)  = (X,!  xX,2  )-X,3.  Note  the  special  splitting  of  the  volume  inte- 
gral in  two  parts.  Through  this  split  it  is  possible  to  apply  the  same  interpo- 
lation functions  with  respect  to  the  reference  surface  for  the  multi-director 
formulation  and  for  the  one-director  kinematics.  Linearization  yields 

DG{ u,  77)  • Au  = f [ f ( 5Ek,CklmnAEmn  + SklA8Eki ) Jd&3]  dQ . nq) 
(ft)  (e*) 
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5.  Finite  Element  Discretization 

The  four  node  quadrilateral  Ql-element  with  isoparametric  bilinear  shape 
functions  is  used  in  the  reference  surface  for  all  kinematic  quantities.  For 
one-director  kinematic  shell  elements  with  5 or  6 d.o.f.  per  node  and  for 
multi-director  kinematics  shell  elements  with  more  than  6 d.o.f.  per  node 
are  applied.  For  coupling  both  types  of  finite  elements  a special  transition 
element  is  used. 

5.1.  ONE-DIRECTOR  ELEMENT 

The  approximation  of  geometry  and  displacements  reads 

X0a  = E Nk{€>  V)  X0A' » and  v„  = E NK{Z,rj)\K,  (20) 

K= 1 K=  1 

where  vk  = [u0;  A/3]£-  is  the  nodal  displacement  vector;  it  consists  of  the 
displacement  components  of  the  reference  surface  and  the  increment  of  the 
rotational  vector.  Depending  on  the  place  of  the  node  (within  flat  areas 
or  at  intersections)  the  rotational  vector  is  parametrisized  with  respect  to 
the  local  basis  (5  d.o.f.)  or  the  global  basis  (6  d.o.f.).  Following  eqs.(l)-(3), 
all  necessary  kinematic  values  are  computed.  With  the  differential  operator 
matrix  B,  see  [5],  we  get 

[ 5e  , 6k  , <$7  ]T  = E ■ (21) 

K=  1 

Applying  the  virtual  work  principle  (16)  and  linearization  (17)  yields 
G(v,  <5v)  = E SvK  gK  » with  gk=  I (B^S  - NKt)dCl , 

K~1  4 4 

DG(v,<5v)Av  = E E <5v£Kja,AvL  , (22) 

K=l  L=l 

with  K kl  — f (Bj-CB l + G Ki)d£l  ■ 


5.2.  MULTIDIRECTOR  ELEMENT 

One  major  advantage  of  the  multi  director  formulation  in  contrast  to  a 
standard  3D-finite  element  is  the  2D-data  structure  with  bilinear  shape 
functions  with  respect  to  the  reference  surface,  which  allows  for  an  easy 
coupling  with  standard  shell  elements.  The  approximation  of  the  geometry 
and  the  displacement  field  of  layer  i,  eqs.  (7)-(5),  follows  from 

xh  = E v)  Xk  , Xka  = X0jc  -f  Q3 13 k , 

K=  i ^23) 

< = E nk(Z,v)&uk. 

K= 1 
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The  number  of  components  for  ur-  is  3 (N  + 1)  -f  3N(m  — 2).  Again,  with 

4 

the  differential  operator  B-matrix,  see  [4],  <$E  = J2  B k £ur-.  Putting  into 

K= 1 


(18)  and  (19)  yields 

4 

G{u,rj)  = SnKGK, 


K=i 

with  Gr-  = f £ f Bl-SJ^dCj]dQe-  f NKpJ dQ.ea , 
(Oe)  j=1(©3)  (A«) 

net  net 


(24) 


£>G(u,77)Au=  SuK  K*l  Aiir  , 


K=1L=1 

withK kl=  f [ 

(0.)  j=1(©3) 


(B^PBi  + G^J^dOe 


5.3.  TRANSITION  ELEMENT 

For  coupling  standard  one-director  with  multi-dirctor  shells  a transition 
element  is  used.  To  prevent  sensible  disturbances  of  the  3D-stress  state 
each  layer  i has  to  be  allowed  for  a constant  thickness  strains  Eh  at  the 
multi-director  side  of  the  transition  element  and  C°-continuity  only  in  the 
mid-surface  at  the  other  one-director  side,  see  [6].  Such  a thickness  jump 
of  the  deformed  laminate  is  admitted  there. 


5.4.  SHEAR-LOCKING 

All  presented  elements  use  a special  interpolation  of  transverse  shear  strains 
to  prevent  “shear-locking”,  namely  an  assumed  natural  strain  (ANS)  method, 
see  [3], 


' 74  ' 

1 

' (1  - v)l 4B  + (1  + vh ZD  ' 

. lv  . 

~ 2 

. (1  - OlvA  + (1  + fbrjC  . 

In  this  formulation  shear  strains  in  the  mid-side  nodes  [M  = A,  B,  C , D) 
are  computed  by  the  standard  bilinear  shape  functions  Nk(€,v)  °f  the 
element. 


6.  Example 

Steel- Cord-Reinforced  Rubber  Beam  under  vertical  load. 

Cross-sec.:  width/height  = 100/20;  Layersequenz  [iso/20°/iso/ - 20°/iso]. 
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Equat.  CPU 
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1.  One-Dir.-Shell 

2.  Multi-Dir.-Shell,  N=15 


3.  Coupling  1.  and  2. 

Figure  3.  Composite-Beam  in  A and  B 


Figure  4-  Displacement  curves 
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